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This paper shows a simple implementation of the Histogram Method for extrapolations in Monte 
Carlo simulations, using the moments of the operators that define the energy, instead of their 
histogram. This implementation is suitable for extrapolation over several operators, a type of 
calculation that is hindered by computer memory limitations. Examples of this approach are given 
for the 2-D Ising model. 
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I. INTRODUCTION 

The calculation of the expectation value of some operator $ in a system in contact with a heath-bath of temperature 

T 

3 i ( )T " E>eM-E t /k B T) ' (1) 

is the usual task of equilibrium statistical mechanics. Here E t and $j are, respectively, the energy and the value of 
<j> associated to the z-th configuration, and the sum runs over all possible configurations. These expectation values 
are usually impossible to calculate analytically, but there are many techniques to approximate expectation values 
in the canonical ensemble. Some of them are quite powerful and general, like high-temperature expansions and 
£j ■ renormalization group techniques. 
— i \ A different approach, quite successful in many applications, is the numerical simulation ("Monte Carlo Simulation") 
of the ensemble pIUlI- This paper is concerned with the Histogram Method, introduced years ago to extrapolate the 
results of a Monte Carlo simulation conducted at a point in parameter space — for instance, some given temperature 
and magnetic field — , to a range of those parameters. It was formulated by by Ferrenberg and Swendsen 0, 0> 
although there are some earlier proposals 6]. For completeness, a short description of this technique follows: Assume 
for instance that one has a model where the only control parameter is the temperature, and wants to compute the 
behavior of some function / of the energy E. At some given T one performs a Monte Carlo simulation using any given 
algorithm that generates Q configurations with the correct probability, given by their Boltzmann weight, measures 
f(E) for each configuration, and calculates finally the average value (f(E)). This average may be taken directly, or 
going through the preliminary construction of a normalized histogram Wt(E) for the energies found in the simulation. 
One can write then 

1 Q 

</(£)>t « n E K E J = E W T (E)f(E), (2) 
which is approximation to Eq. JJ) with <£> = f(E). This equation can also written in terms of the density of states 

O g(E) 

, f( ^ _ J2E9(E)f(E)eM-E/(k B T)) 

{f{E))T ~ J2 E 9(E)eM-E/(k B T)) • (3) 

Comparing Eqs. J2J and © it is clear that Wt(E) is proportional to g(E) exp(— E / (k B T)), and so an approximation 
for the density of states — ignoring normalization for the moment — is given by g{E) « Wt{E) exp(E/(k B T)). Notice 
that g(E) is a property of the system itself, independent of temperature. 

With an approximated density of states at hand the results of the simulation are extended to other temperatures. 
Eq. JSJ written for a different temperature T' gives 

_ J2 E g(E)f(E)eM~E/(k B T)) 



(f(E)h 



Y jE g{E)eM-E/{k B T')) 

(4) 



J2e W T (E)f(E) exp(E/(k B T) - E/(k B T')) 



J2 E W t {E) exp(E/{k B T) - E/(k B T')) 
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and one says that the histogram Wt(E) has been reweighted. 



II. EXPANSION IN MOMENTS 



In general, with the Histogram Method results obtained at a given point in the space of control parameters can be 
extrapolated to a neighborhood of that point. However, the implementation of the method can be cumbersome for 
models with more than one parameter, because the size of the histograms generated becomes too large. For instance, 
an Ising model in a simple cubic lattice, with nearest-neighbor, next-nearest-neighbor and magnetic couplings has 
a 3D parameter space. When applying the histogram method to this system one faces the problem of handling a 
histogram with a total number of bins that rapidly becomes impossible to accommodate in the computer's physical 
memory. For instance, for a very modest size L — 16, one needs to store some 16 3 «4x 10 3 possible values for each 
of the 3 couplings, and this gives a histogram with around 64 x 10 9 bins. Working with 4 Byte integers (which may 
be insufficient for very large runs) gives a memory requirement around 250 GBytes, realizable in large workstations 
and supercomputers, but usually beyond the reach of the commodity machines used by many scientist today Q. 

Among possible solutions one could coarse-grain the binning of the histogram (this is unavoidable if one is dealing 
with continuous variables, like in the XY or Heisenberg Models); one could also store a record of the operators of 
interest for all configurations generated by the Monte Carlo algorithm (this has been recommended for Hamiltonians 
with continuous variables |2|, |8()- Here I propose a different option that arises from two simple facts: one, most 
of the quantities one is interested in calculating can be expressed in terms of moments of the same operators that 
constitute the Hamiltonian. Two, the reweighting function can be easily expanded in a power series. Then, the 
physical quantities one wants to estimate at different values of the couplings can be obtained as power series in the 
increments of said couplings, using the moments of the operators obtained at some fixed point in parameter space 
as coefficients. A similar approach was proposed in Ref. using cumulants instead of moments; This approach via 
moments is bit simpler. Here I also show some details of implementation needed to insure larger ranges of numerical 
convergence. 

To fix ideas, I shall use a simple example, namely, the Ising Model with zero magnetic field. The Hamiltonian is 

H = - J ]T S a S b . (5) 

<a,b> 

Assume now that one calculates the expectation value of some power of the magnetization, say (A/ 2 ), at some 
given temperature T. Define for convenience the adimensional energy operator 6 = E< a b> SaSb and the couplings 
K = J/k B T. From Eq. (QJ one gets 

E, exp(^9 l ) 

Now, if one wants to calculate this same expectation value at some other temperature, one could store the histogram 
for and then reweight it. But one can also write the desired expectation value 

, M 2, = E t Mfexp(iT9 z ) = Mj exp(if 6Q exp((iT - K)Qj) 

Eiffltp^eo E i ex P (^e I )exp((iv'-A')e l ) ' {> 

and expand the exponential containing K' — K in a Taylor series 

/wax _ Ei M? exp(ge0 EZo (j? ~ K) n Q?/nl 

{ )K ' ^MK^ZZoiK'-Ky^/n! ' W 
For finite systems the order of the sums can be safely interchanged, giving 

iM 2, = EZo [Cg - k'TIA mm ex P (if e.) 

[ }K ' E,T=o P" - K ') n M E, ©r ex P (A^) • 1 > 

One can now divide numerator and denominator by the partition function Z(K) and find that the terms at the end of 
both expressions are the expectation values (moments, in short), of Af 2 0™ and <d n , respectively, taken at K. Denoting 
8K = K — K 1 one gets 

{M )ki = ^mwfcM ■ (10) 
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Since the expansion of the exponential function is convergent everywhere, there is no a priori limitation for this 
approach. But it is clear that one can never actually calculate all the required moments, and a truncation needs to 
be applied. This immediately introduces very strong bounds in the applicability of the method, since for large 5K 
the two sums in the previous expression require the inclusion of higher and higher moments if reasonable results are 
going to be obtained. A discussion about the number of moments needed to insure certain range of convergence is 
given later on. 

Once the series are truncated, one faces a more serious difficulty: consider again the M 2 example, and assume that 
a simulation has been conducted in a lattice with TV spins. Since M is an extensive quantity, (M 2 ) is of order TV 2 , and 
the range of K where some degree of numerical convergence can be achieved is quite small. To find a way around this 
problem — at least partially — , begin by working with densities instead of with the original operators. Going back to 
our example, define m = M/N and 9 = Q/N. The previous expression is rewritten as 

_ EZo(N6Kr(m 2 e") K /n! 

Now it is clear that, since (9) is of order one, one may expect numerical convergence of the strongly truncated series 
up to 

Aff«i - AT«^, (12) 

which is a very narrow range; here 'strongly truncated' means a series with a few terms, say, less than 20. This is 
actually a very conservative bound, since the n! in the denominators improve numerical convergence (and guarantees 
it for the infinite series). 

The convergence of the truncated series can be improved if one notices that the desired moments are taken from 
a sharply peaked distribution. The combination of the fast increasing density of states <?(0) and the fast decreasing 
Boltzmann weight implies the existence of a narrow distribution Wk(9) — that is, the histogram in 9 — , centered in 
some value 9 C . Going back now to Eq. ifSjl. one may see that the zero values of the operators appearing in the energy 
(and therefore in the energy density) can be easily shifted. For instance, if one shifts the density 9 by some value 9r, 
the expectation value for m 2 becomes 

2 j: i m 2 exp(NK9 i ) _ m 2 exp{NK{9 t - 9 R )) 
Km ' K E I exp(™ i ) ^expiNKiOi-en)) ' 1 ' 

since the factors of exp(—NK9 R ) cancel in the fraction. Now, if the value one chooses for 9r is close to the center 
of the distribution, the moments one gets for 89 = 9 — 9 r are going to be numerically small — and become smaller 
the higher the moment and the larger the lattice — , simply because of the narrowness of the originating distribution. 
Retracing the steps taken going from Eqs. Q to (|ll)fl one gets the final expression for (m 2 )K' as 

/ m 2v _ Zt (N6Kr(m 2 (69) n )K/n! 

where both numerator and denominator series have been truncated to n* terms. Now the fast growth in the NSK 
coefficients is partially balanced by a fast decrease in the values of ((S9) n )K and (m 2 (69) n ) k , and in this way the 
applicability of the method is extended to a much wider range in couplings. It is not too difficult to calculate the 
range of expected convergence of the extrapolation. In fact, for temperatures away from the critical point the width 
of Wk iff) is order 1/y/N, from where one gets a range of convergence for the extrapolation of 

1 T 2 
AK w — == -> AT w — =, (15) 

which comes from estimating N AK A9 1. Close to the critical point the width of Wt(9) increases, and scales as 
A9 w N a / { 2vd ) j y/N. From here one estimates then a range of convergence 

T 2 

AT m - . (16) 

Here a and v are the exponents for specific heath and correlation distance close to the critical point, and d is the 
dimensionality. For a small a this is still a much wider range than the one gotten when no shifts are used. Notice 
that for the 2D Ising Model one replaces N a ^ vd ^ by In TV and gets an scaling 

T 2 

ATk ■ =. (17) 
V/VlniV 
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Notice finally that the example given here, that extrapolates for m 2 , can be easily extended to any other lattice 
operator. 



A numerical test of the method was done for the two-dimensional Ising Model in square lattices of sizes 16, 32 and 64, 
with the Hamiltonian given in Eq. (JSJ . Two Monte Carlo simulations were carried out using the Wolff algorithm j 
at a temperature of fc^T '/ 'J = 2.27, close to the critical temperature for the model. All moments of the form (m l (80) n ) 
and ((S\m\) 1 (56) n ) for < I, n < 16 were recorded. First a short run with no shifts in either 6 or \m\ was done, using 
2 x 10 6 , 3 x 10 6 and 4.4 x 10 6 Wolff iterations for L = 16, 32 and 64, respectively, after discarding transients of 120, 
170 and 270 iterations. The running time for this set of simulations was about 1,800 seconds in a 2.4 GHz Pentium 
4 processor. The exact solutions for the energy and the specific heath were obtained from the analytic free energy 
for finite lattices found by Kaufman Fig. (1) shows the results of the extrapolation for the adimensional energy 
density (9) using either 14 or 15 moments |12j, when no shift in 9 has been implemented, compared with the exact 
results. It is clear that the range of applicability of the extrapolation is extremely narrow, as expected from Eq. I|12|l. 
That estimate gives here AT 32 = 0.00 50 and AT 64 = 0.00 1 3. This estimate is quite conservative, and the figure gives 
ranges of numerical convergence which are 2 to 3 times longer. 

The second simulation was a larger run where the results for (9) and (\m\) from the first were used as reference 
values 9r and \m\n for the shifts. The total numbers of Wolff iterations were 4 x 10 7 , 6 x 10 7 and 8.8 x 10 7 for 
L = 16, 32 and 64, respectively, after transients of 120, 170 and 250 iterations were discarded. The total time 
for this set of simulations was about 37, 000 seconds with the same processor. Data were divided in 20 blocks in 
order to generate error estimates. Extrapolation for energy and specific heath were calculated. The results for the 
energies are given in Fig. (2), those for specific heath in Fig. (3), and in both figures a comparison with the exact 
results is given. It is clear that the range of applicability of the extrapolation is larger, and actually a bit larger 
than the estimation made in Eq. 1)170. which gives AT32 = 0.061 and AT^ = 0.028. Two important points should 
be remarked: First, for (9) the direction in which the extrapolated curves deviate from the correct results, for low 
T, depends on the number of moments taken into account; they deviate upwards when one uses 14 moments, and 
downwards when using 15 moments. A similar behavior is found for c, except that now the two extrapolations deviate 
in different directions at both ends of the range of convergence. This gives a very simple way of bounding the range 
of convergence of the algorithm. Second, the statistical errors in the extrapolations grow as one moves away from 
the simulated temperature, but the effect of this growth is smaller than the effect of the change in the number of 
included moments. The behavior of both types of errors are given in Fig. (4), which shows the difference between 
extrapolated and exact energies for L = 32. The behavior for errors is similar for the specific heath. Notice that the 
statistical errors in the extrapolated results actually become smaller for temperatures a bit below the point where the 
simulation was carried out; this curious phenomenon has been studied in full histogram extrapolations |8(- 

The moments of \m\ were used to generate an extrapolated estimation of the susceptibility x', defined by x' = 
NK'((\m\ 2 ) — (|m|) 2 ) = NK'((m 2 ) — (\m\) 2 ). This expression is used instead of the true su scep tibility x which in 
a numerical experiment does not give the expected peak close to the critical temperature [13}. Fig. (5) shows a 
comparison between the extrapolated x' an d values obtained in other individual simulations. These were obtained 
at their nominal temperatures using the Wolff algorithm. It is remarkable how the extrapolated results manage to 
reproduce reasonable well the peak in the susceptibility. Here one can also notice that the order of the approximation 
again decides in which direction the extrapolated results deviate form the actual values, and so a simple comparison 
between the 15- and 16-moments expansions gives bounds for the region of convergence of the method. 



As mentioned before, this algorithm becomes attractive especially in cases where one has to deal with Hamiltonians 
composed of several operators, where the sizes of the histograms needed for reweighting may overflow the available 
memory. As a simple example consider again the 2-dimensional nearest-neighbor Ising Model, but now with a magnetic 
field. The Hamiltonian is 



A. Example: the 2-D Ising Model: extrapolations at h = 



B. Example: the 2-D Ising Model with magnetic held 




(18) 



<a,b> 



a 



giving a Boltzmann weight 



exp(-H/k B T)) = exp(KO + hM) 



(19) 
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where the additional definitions M = J^a^a and h = H^/UbT have been introduced. Consider now a simulation 
carried out at some temperature and magnetic field. Assuming that one wants to extrapolate the expectation value 
of some operator $, one gets 

= ^ eMKQ> + hM t ) eM(K' - K)Sj + {h' - h)Mj) 
E,rx,,(A(-), • //.\/,;cx,.((A'' A;(-), • (h> h)Mr 



As before, it is better to change all operators into their densities, defining to = M/N and <j> = $>/N. Also, the density 
6 should be shifted so that its distribution is centered close to zero (this is unnecessary for to, since its distribution 
is symmetric around to = 0). Performing these operations and the Taylor expansions for the exponential containing 
SK and Sh, one gets, repeating the same steps that gave Eq. I|14|) . the following expression 

Y.r=vY.n=*{N5K) l (N8Kr{d>{59) l m") K , h /{l\n\) 
[<P)K ' H ' EZoEZ (N6Ky(N6hr((69y m -) KM /(Unl) ' 1 > 

One should pay attention here to the fact that for low temperatures the histogram Wk,Ii(6, m) becomes bimodal in to, 
and the assumption of a distribution with a single narrow peak is no longer valid. However, close to the critical point 
the width in m of such histogram remains small, and one can still get by with the first few moments. Occasionally, 
it may be convenient to work in terms of |m|, whose distribution remains unimodal. 

The data obtained from the previous simulations were now used to generate the behavior of the magnetization and 
the susceptibility at non-zero values of H . For the magnetization one gets, after truncation 

(m) = EloE:'= (NSKY(Nhr((89y m ^) K j l=0 /(lM) ^ 
El EnU N ^y(Nh) n (m i mn) K<h=0 /(l\rd) ' 

an analogous expression is obtained for (m 2 ) k> ,h' > anc ^ from here the true susceptibility can be computed. The results 
are shown in Fig. (6), which shows: (a) magnetization vs. H for L = 16, 32 and 64 at the critical temperature 
T c = 2/ln(l + V2) = 2.269185..., and (b) magnetization vs. H for L = 32 at T = 2.20, 2.27 and 2.34. In all cases 
the magnetization has been extrapolated from the simulations that were carried out at H — and T = 2.27 '. In a 
slight departure from what was done before, here the denominator in Eq. (|22|l was calculated using 16 moments while 
the sums in the numerator were truncated to 15 or 14 moments. The individual points have been calculated using a 
modified version of the Wolff algorithm, were the acceptance ratio depends on the change of energy due to a cluster 
flipping in the presence a magnetic field H. It is clear from the figure that the expansion in moments reproduces quite 
well the behavior of the magnetization for each temperature and lattice size, and in particular it manages to show 
the large growth of to with H as T is reduced. The splits in the extrapolation curves correspond to the separation of 
the 15- and 14-moment extrapolations, and mark the end of the ranges of convergence. For small lattices these splits 
do not appear in the H range tested here. 

Finally, Figs. (7) and (8) show the true susceptibility x as a function of H , for T — T c and L = 16, 32 and 64, and for 
L = 32 and T — 2.20, 2.27 and 2.34. It should be noticed that the extrapolation manages to cover quite well the whole 
peak in \- The results for L — 32 and different temperatures also show an excellent agreement between extrapolations 
and individual simulations, and show how the method can really extrapolate in more than one parameter. Notice 
that, as expected, \ grows as T is reduced, even below T c ; for H = the true susceptibility is just proportional to 
(rn 2 ). Otherwise, the behavior the extrapolated quantities vis a vis the individual simulations at H ^ is in all 
respects analogous the results found before: a very good correspondence for small H — that is, close to where the 
simulation was carried out — , followed by a large deviation, which depends on the number of moments included in 
the extrapolation. 



III. CONCLUSIONS 



This paper shows how to implement the histogram method for extrapolation of results of a Monte Carlo simulation 
using the moments of the histogram. This approach has several advantages over the direct method — histogram 
construction and posterior rcwcighting — , and over the method of storing configurations for their reweighting (named 
"histogram on the fly" in Ref. |8|). To start with, the resulting expressions for the extrapolated quantities are given 
by very simple and conceptually appealing formulas. Second, the ranges of applicability of the method become 
evident simply by changing the number of moments included in the extrapolation. Third, the amount of computer 
memory and physical storage needed are so small that one may without any problem generate several repetitions of 
the simulations so as to generate in a simple way the error estimates for the extrapolated quantities. And finally, this 
approach eliminates completely the need to choose a binning size in cases of continuous variables. 
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One needs however to balance these benefits against the cost of the extra approximation involved in the method. 
After all, replacing a full histogram for its first few moments necessarily reduces precision. How many moments 
one really needs to keep in any given simulation so that not too much information is lost is an issue that has to be 
considered with some care. On the one hand, it is clear that increasing too much the number of stored moments not 
only defeats one of the motivations for this approach, which is to work with a limited memory, but also necessarily 
runs into the limits of reliability imposed by the statistical errors of the simulation. On the other, not keeping enough 
moments implies a waste of simulation time. As a first approach to the answer to this question one can consider 
the following estimation, done here for a one-coupling Hamiltonian: consider the preliminary run needed for the 
estimation of 9 C . This preliminary run can also be used to obtain a rough estimate of the width <jg of the distribution. 
Now, the order of magnitude of the moments ((S6) n ) will be around <jg, and so the terms needed in an estimation 
with a shifted coupling look like (see Eq. i|14|)t 



It is not difficult to show, using a steepest descent approximation, that as a function of n this expression behaves as 
a Gaussian centered in n max = AK tj$, with a width given by ^/n max and height exp(n max ). Therefore one gets 



The conclusion is then the following: given an initial estimation of the width ag of the distribution, and assuming 
that a certain maximum extrapolation range AK is going to be used, the main contribution to the extrapolation 
comes from the moments with n around n max = AKag. Besides, one finds that the moments with n such that 
n — "-max 3> ^Timax are basically irrelevant. The numerical results shown here display a much larger dependence on 
the number of moments used than on the statistical spread of the data, suggesting that several more moments may 
have been used in the extrapolation before reaching the limits given by statistical spread. 
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FIG. 1: Energy per site vs. Temperature for the 2-D Ising Model, for L = 32 and L = 64. The dotted line gives the exact results 
obtained from the Kauffman solution. The solid lines give 14- and 15-moment extrapolations using values from a simulation 
at T = 2.27 (indicated by the vertical line). Here the moments are taken without using any shifts [Or — 0). 
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FIG. 2: Energy per site vs. Temperature for the 2-D Ising Model, for L = 32 and L = 32. The dotted line gives the exact 
results obtained from the Kaufman solution. The thick solid lines give 14- and 15-moment extrapolations using values from a 
simulation at T = 2.27 (indicated by the vertical line). The thin solid lines give the error bars for the extrapolations. Here the 
moments are taken using values of 9r of —1.432094 for L = 32 and —1.422108 for L — 64. 
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FIG. 3: Specific heath per site vs. Temperature for the 2-D Ising Modei, for L — 16, 32 and 64. The dotted lines give the 
exact results obtained from the Kaufman solution. The thick solid lines give 13- and 14-moment extrapolations using values 
from a simulation at T = 2.27 (indicated by the vertical line). The thin solid lines give the error bars for the extrapolations. 
Here the moments are taken using values of 8 R of -1.451724 for L = 16, -1.432094 for L = 32, and -1.422108 for L = 64. 
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FIG. 4: Error in the extrapolation of the energy density for L = 32. The solid line gives the difference between the extrapolated 
energy and the exact Kaufman solution. The error bars correspond to the extrapolation. Here 15 moments were used. Here 
R = -1.432094. 




FIG. 5: Susceptibility calculated using \m\, for the 2-D Ising model with L — 16, 32 and 64. The points are from 
independent simulations, carried out at their nominal temperatures, an the dotted line is given only as a guide to the eye. The 
solid lines correspond to extrapolations from a simulation at T = 2.27, using 15 and 16 moments. 
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FIG. 6: Magnetization vs. H for the 2-D Ising model. The lines give the values extrapolated from a simulation at T = 2.27, 
H — 0, and the points, with their error bars, correspond to individual simulations carried out at their nominal T and H values, 
(a): Results for T — T c = 2.269185 . . ., with sizes L — 16 (triangles), 32 (squares) and 64 (round marks). The split visible for 
the L = 64 line correspond to the separation between the 14- and the 15-moments extrapolations, (b) Results for L — 32, with 
temperatures T = 2.34 (triangles), T = 2.27 (squares) and T — 2.20 (round marks). In both cases the extrapolated statistical 
errors come out smaller than the thickness of the lines. 
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FIG. 7: True susceptibility \ vs. H for the 2-D Ising model. The lines give the values extrapolated from a simulation at 
T = 2.27, H — 0, and the points, with their error bars, correspond to individual simulations carried out at their nominal T and 
H values. Results shown here are for T — T c — 2.269185 . . ., with sizes L — 16 (triangles), 32 (squares) and 64 (round marks). 
The split visible for the L — 64 line correspond to the separation between the 14- and the 15-moments extrapolations. The 
extrapolated statistical errors come out smaller than the thickness of the lines. 
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FIG. 8: True susceptibility \ vs. H for the 2-D Ising model. The lines give the values extrapolated from a simulation at 
T — 2.27, H — 0, and the points, with their error bars, correspond to individual simulations carried out at their nominal T 
and H values. Results shown here are for L = 32, with temperatures T = 2.34 (triangles), T = 2.27 (squares) and T = 2.20 
(round marks). The extrapolated statistical errors come out smaller than the thickness of the lines. 



